Friedel Oscillations疑问

问题

动量空间中电势为

其中

所以

\begin{equation} \label{eq:1} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 {\varepsilon(\vec{q})}} \end{equation}

Thomas-Fermi 近似的结果

Thomas-Fermi 近似给出

将上式代入 $\varphi(\vec{q})$ 有

\begin{equation} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 \left(1 +\frac{q_{TF}^2}{q^2}\right)} =\frac{-e}{\varepsilon_0 V \left(q^2 +q_{TF}^2\right)} \end{equation}

对其 Fourier Transform

对

在 $q \rightarrow +\infty$ 展开(与在 $q_{TF}\rightarrow 0$ 时对 $q_{TF}$ 展开相同):

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
s = sym.series(q/(q**2+q_tf**2),q,sym.oo,5)
print(sym.latex(s))

- \frac{q_{TF}^{2}}{q^{3}} + \frac{1}{q} + O\left(\frac{1}{q^{5}}; q\rightarrow \infty\right)

计算 leading order

这正好是精确结果

在 $q_{TF}\rightarrow 0$ 时对 $q_{TF}$ 展开的结果的 leading order 相同.

问题1

Thomas-Fermi 近似结果与精确结果的 $q_{TF}^2$ 的系数不同?

答: 因为展开后的结果中出现了发散项, 所以问题本身就不能展开.

问题2

为什么要对 $q$ 在 $q\rightarrow +\infty$ 展开? 积分的区间不是整个实轴吗?

答: 我之前对展开的目的理解有误. 展开的目的不是为了一阶一阶地计算积分. 展开目的是为了数值的计算积分. Thomas-Fermi 的积分本身可以精确解出, 没必要数值积分, 因此不必做展开. 展开的意义会在下面 RPA 的计算过程中体现.

RPA 的结果

RPA 给出

将上式代入 $\phi (\vec{q})$ 有

\begin{equation} \varphi (\vec{q})= \frac{-e}{\varepsilon_0 V q^2 \left(1 +\frac{q_{TF}^2}{q^2}\right)} =\frac{-e}{\varepsilon_0 V \left(q^2 +q_{TF}^2 g\left( \frac{q}{2k_F} \right) \right)} \end{equation}

其中

对 $\phi (\vec{q})$ 作 Fourier Transform

级数展开

对

在 $q\rightarrow +\infty$ 展开

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
kf = sym.Symbol('k_F')
u = q/(2*kf)
*g=1
g = sym.Rational(1,2)*( 1+(1-u**2)/(2*u)*sym.log((u+1)/(u-1) ) )
s = sym.series(q/(q**2+q_tf**2)*g,q,sym.oo,10)
print(sym.latex(s))

- \frac{\frac{256 k_{F}^{8}}{63} - \frac{64 k_{F}^{6} q_{TF}^{2}}{35} + \frac{16 k_{F}^{4} q_{TF}^{4}}{15}}{q^{9}} + \frac{\frac{64 k_{F}^{6}}{35} - \frac{16 k_{F}^{4} q_{TF}^{2}}{15} + \frac{4 k_{F}^{2} q_{TF}^{4}}{3}}{q^{7}} + \frac{\frac{16 k_{F}^{4}}{15} - \frac{4 k_{F}^{2} q_{TF}^{2}}{3}}{q^{5}} + \frac{4 k_{F}^{2}}{3 q^{3}} + O\left(\frac{1}{q^{10}}; q\rightarrow \infty\right)

问题3: 接下来该怎么做? 所有的展开项代入积分都是发散的.

答: 展开算错了. 程序中把公式抄错了. 看来最好还是手算 check 一下比较好, 不然不会出现这么低级的错误. 级数展开的计算能力还有待提高...

正确的级数展开

import sympy as sym
q = sym.Symbol('q')
q_tf = sym.Symbol('q_TF')
kf = sym.Symbol('k_F')
u = q/(2*kf)
g = sym.Rational(1,2)*( 1+(1-u**2)/(2*u)*sym.log(sym.Abs((u+1)/(u-1)) ) )
s = sym.series(q/(q**2+q_tf**2*g),q,sym.oo,10)
print(sym.latex(s))

- \frac{1}{q} - \frac{4 k_{F}^{2} q_{TF}^{2}}{3 q^{5}} - \frac{16 k_{F}^{4} q_{TF}^{2}}{15 q^{7}} - \frac{64 k_{F}^{6} q_{TF}^{2}}{35 q^{9}} + O\left(\frac{1}{q^{10}}; q\rightarrow \infty\right)

数值计算 RPA 的结果

问题分析

首先, 积分

难以解析地计算, 因此想要数值地计算.

其次, 被积函数收敛不够快.

经过前面对被积函数级数展开的分析, 发现被积分函数在远处的 leading order 为

\begin{equation} \label{eq:leading} \frac{\sin(qr)}{q} \end{equation}

也就是说, 被积函数在远处的贡献主要来自于 leading order 项 (\ref{eq:leading}), 而且这一项是可以解析地计算出结果的. 如果我们从被函数中把这项减去, 那么被积函数在远处就近似为 $0$ 了, 也就是说积分会更快地收敛. 这就解决了积分收敛不够快的问题.

构造新的被积函数

构造函数

\begin{equation} F(q) = \sin(qr)\cdot \frac{q}{q^2+q_{TF}^2g( \frac{q}{2k_F} )} - \frac{\sin(qr)}{q} \end{equation}

新的 $F(q)$ 函数会很快地收敛.

原来的积分就可以分解为

数值积分

下面数值地计算

python 程序如下

#导入数值计算包, 画图包, 积分包
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

#设置物理参量
qtf = 3/4
kf = 1/2

#数值积分的区间:(0,end)
end = 5

#定义函数 phi_r, g(u), phi_q, lo(q)是phi_q 的 leading order, F(q)是从 phi(q) 里边减去 leading order
def phi_r(r):
    def g(u):
        return 1/2*(1+1/(2*u)*(1-u**2)*np.log(np.abs((1+u)/(1-u))))
    def phi(q):
        phi = np.sin(q*r)*q/(q**2+qtf**2*g(q/(2*kf)))
        return phi
    def lo(q):
        lo = np.sin(q*r)/q
        return lo
    def F(q):
        F = phi(q) - lo(q)
        return F
    (Fres,err) = integrate.quad(F,0,end)
    phir = Fres + np.pi/2
    return phir


#计算N个不同的 r 的取值的数值积分 
N = 100


#画六个图中的第一个图
plt.subplot(321) #(321)的意思是整个图是3行2列,这个图画在第1个位置
a = 0.1 #画图的区间为[a,b],这里第一个图为 r 取 [0,10]
b = a+10
r = np.linspace(a,b,N)

#进行 N 次数值积分, 然后把 N 个点的数值积分的结果保存在数组 s 中
s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
#画图, 指定线的颜色和图例
plt.plot(r, s, color="red", label="Numerical Result")
#我们想要将数值的结果和 Friedel Oscillations 的结果做比较, 但它们之间会差一个常数k, 现在把这个常数 k 取 N 次数值积分中第一个点的比值
k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r) #这是 Friedel Oscillations 的结果
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation") #画出 Friedel Oscillations 的结果

plt.xlabel("$r$") #设置横坐标的物理量
plt.ylabel("$\phi(r)$") #设置纵坐标的物理量

#以下三行是为了让坐标的数值以科学计数法显示
ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。
plt.legend()



#与画第一个图同样的方法, 画出剩下的5个图, 它们的区别在于 r 的区间不同
plt.subplot(322)
a = 10
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。

plt.subplot(323)
a = 20
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(324)
a = 50
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(325)
a = 100
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。


plt.subplot(326)
a = 200
b = a+10
r = np.linspace(a,b,N)

s = np.linspace(0,0,N)
for i in range(N):
    s[i-1] = phi_r(r[i-1])
plt.plot(r, s, color="red", label="Numerical Result")

k = phi_r(a)/(1/(a**3)*np.cos(2*kf*a))
fo = k*1/(r**3)*np.cos(2*kf*r)
plt.plot(r,fo, "--", color="blue", label="Friedel Oscillation")

plt.xlabel("$r$")
plt.ylabel("$\phi(r)$")

ax = plt.gca()  # 获取当前图像的坐标轴信息
ax.yaxis.get_major_formatter().set_powerlimits((0,1)) # 将坐标轴的base number设置为一位。



plt.suptitle("The Results of RPA") #画出整个图的标题
plt.show() #显示图片

#+RESULTS: : None

结果对比

将数值结果与 Friedel Oscillations 的结果

进行对比

可以看出, Friedel Oscillations 的结果在 $r$ 较大时的结果比较好.

参考文献

Friedel Oscillation 的原始文献 the shielding of a fixed charge in a high-density electron gas http://www.doc88.com/p-9512851691956.html

Wolfgang Nolting, Fundamentals of Many-body Physics

致谢

感谢导师 Ran Qi 的指导

感谢 Fan Yang 师兄的讨论